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Abstract 


Interpolation at grid boundaries is studied for the purpose of solving 
partial differential equations using either implicit or conservative explicit 
finite-difference methods on multi -component overlapping grid systems. 
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Introduction A multi-component grid system, in which several computation 

grids, are used, is required in the numerical solution of many fluid dynamics 
problems involving flow within or about a complicated geometric configuration. 
From a grid construction point of view, the simplest procedure is to generate 
each component grid independently with a sufficient overlap so that information 
can be transmitted from one grid to the other. The development and analysis 
of solution procedures on this type of grid system was studied by Starius 
[8,9]. The practical application of the method to the solution of problems in 
computational fluid dynamics was demonstrated in the papers by Atta [1], Atta 
and Vadyak [2], and Thompson [11]. This was followed by further studies on 
interpolation techniques by Kreiss [5] and Mastin and McConnaughey [6]. Each 
successful application, such as the recent results of Steger and Buning [10] 
and Benek, et al [ 3 ], serves to reinforce the need for additional work on the 
implementation of numerical methods on overlapping grids. 

Two popular numerical methods in computational fluid dynamics are the 
alternating direction implicit algorithms and the explicit algorithms derived 
for the solution of conservation laws. Both methods have been used in the 
solution of problems on composite grid systems. However, in each case there 
are properties of the numerical solution which are lost when information is 
transmitted between the individual grids. With the implicit algorithms, there 
is no technique for generating advanced solution values at all boundary points 
of each component grid. Lagging some of the boundary point values can lead to 
a loss of accuracy in the solution of transient problems. Of course, the 
temporal step length could be reduced, but that would defeat the purpose of 
choosing an implicit method. It is also possible that lagging may effect the 



stability of the method, although no problems of that kind have been reported. 

A conservative finite-difference scheme is often selected when solving 
partial differential equations in conservation form. Since the classical 
interpolation formulas were not derived with conservation properties in mind, 
their use in finite-difference approximations on composite grids would result 
in the loss of. an exact conservation property. An interpolation scheme for 
conservative finite-difference methods was first proposed by Rai [7]. He 
considered composite grid systems which did not overlap but joined along 
common grid lines. Berger [^] indicated how the method of Rai could be 
generalized and extended to overlapping grids. 

This report will describe ways of eliminating the time lag in implicit 
solutions and will present a general algorithm for constructing conservative 
interface conditions. The variables x and y are used as the spatial variables 
in the partial differential equations. Since, in the general case, one would 
need to solve the equations on a curvilinear grid, all equations would be 
transformed to curvilinear coordinates before applying the finite-difference 
algorithm. In each development the partial differential equation is 
sufficiently generally so that the method can be applied to the original or the 

V . 

transformed equation without change. For simplicity, mixed derivative terms 
which are normally lagged in the single grid case and source terms are not 
included in the partial differential equations. 



Implicit Methods 


The fundamental concepts are quite simple and can be 


demonstrated by considering the one-dimensional equation 

u t = Lu , (1 ) 

where L is a differential operator of the form 

Lu = au x + bu xx . 

Since implicit methods generally require linearization of the difference 
equations, it may as well be assumed that L is linear. Suppose that two grids 
G-| and G 2 are given on the intervals [a,d] and [c,b], respectively, where 

a < c < d < b. 


If M denotes the usual second order difference approximation of L, then the 
Crank-Nicolson equation can be written as 

,n+1 ..n . At /K ,..n+1 

2 


u i = uj + fE (Muf 1 + Mu}). 


( 2 ) 


Here i is the spatial index, n is the temporal index, and At is the step 

length. Now suppose values on G-| and G 2 are known at level n and values at 

level n+1 are to be computed on G 1 . While solution values needed in (2) at 

x=d can be interpolated from G 2 for level n, the corresponding values at level 
\ . 

n+1 are unavailable. If these unknown values are replaced by the values at 
level n, then the local truncation error at the neighboring interior point is 
increased by a term on the order of 0(At /Ax ). The value Ax represents the 
spatial grid spacing on G^ , or the spacing at x=d in the case of a nonuniform 
grid. In any event, when Ax is small, this lagging of solution values will 
seriously degrade the temporal accuracy of the approximation. The error can be 
reduced by following a particular order in updating the solution values at the 
interior grid boundary points. The correct sequence of computations is 
indicated in the following steps. 



1 . 

Calculate 

u n+1 

on 

G 1 

with 

level 

n values at 

x= 

d. 

2. 

Calculate 

u n+1 

on 

G 2 

with 

level 

n+1 

values 

at 

x=c. 

3. 

Calculate 

u n+ 2 

on 

G 2 

with 

level 

n+1 

values 

at 

x=c. 

4. 

Calculate 

u n+2 

on 

Gi 

with 

level 

n+2 

values 

at 

x=d. 


Now the error induced by using the previous value at x=d in step 1 is offset 
by the use of the advanced value in step *4. In fact, the local truncation 
error at the neighboring interior point is increased by a term of order 
OCAt^/Ax 2 ) when the solution is advanced from level n to n+2. The same error 
reduction would also occur at x=c. 

Clearly, this four-step alternating grid scheme is only a partial 
solution. Unless the solution exhibited a linear growth or decay, there would 
still be points with a local truncation error of order one whenever At=Ax. 
However, this does not necessarily mean that the global error in the numerical 
solution would be increased to that order. The actual error in the solution 
would also depend on other factors such as the extent of the overlap. Note 
that the same updating procedure could be applied to implicit methods other 
than the Crank-^Nicholson method, but the reduction in local truncation error 
would not be the same. 

The alternating grid concept has also been used in the development of 
another method for implementing implicit algorithms on composite grid systems. 
This method also alternately employs the forward difference explicit equation 


,n+1 


u? + At Mu? 


(3) 



and the backward difference implicit equation 

u j* + 1 = u^ + At Mu^ +1 . (H) 

The computational sequence is illustrated in the following four-step procedure 
which advances the solution from level n to level n+2. 

1. Calculate u n+1 on G 1 using (3). 

2. Calculate u n+1 on G 2 using (4). 

3. Calculate u n+2 on G 2 using (3). 

4. Calculate u n+2 on G^ using ( -4 ) . 

The method alternates the explicit and implicit calculations in the same 
manner as the well-known hopscotch algorithm. Thus, the name hopscotch will 
be associated with this method. The method has several desireable properties. 
All values needed at the grid boundaries c and d can be computed by 
interpolation from solution values at the correct time level. The overall 
method is second-order accurate in time and unconditionally stable. This fact 
follows by noting that the combined sequence of (3) followed by (4) is 
equivalent to a Crank-Nicolson step with step length of 2At. 

The consequences of lagging solution values at grid boundaries can be 
even more serious for multi-dimensional problems. Suppose for example that 
the operator L in (1) is defined as 

Lu - au x - bu y ♦ cu xx * du yy . 



Let M and M„ denote the difference approximations of the x and y derivative 

a y 

parts of L. If the parabolic equation (1). is solved by an ADI method, such as 
Peaceman-Rachford, the algorithm becomes 


u n+1/2 _ u n + (MxU n + 1/2 + MyU n } f 

u n+1 = u n+1/2 + — (M Y u n+1/2 + M u n+1 ) . 

2 x y 

Now the error that occurs in the first step of the algorithm is further 
magnified in the second step. This argument can be made more precise by 
noting that the Peaceman-Rachford ADI method is a perturbation of the 
two-dimensional Crank-Nicolson method with a perturbation term 

— - M 2 M 2 ( u n - u n+1 ) . 

A lagged value in this term produces a truncation error term on the order of 
0(At^/Ax^r 2 ). The alternating grid procedure would reduce this to 

Jl p p 

0(At /Ax Ay ). The one-dimensional hopscotch algorithm would not be a 
computationally efficient method for two-dimensional problems. However, the 
same effect can be realized by inserting additional steps in the ADI algo- 
rithm. The procedure is again demonstrated using two grids G-, and G 2 . Note 
that equation (5a) can be written as 


(5a) 

(5b) 


,n+1 /2 _ „n ^ At 


u n + — (M x u n + M y u n ) 


u n*l/2 . v n + 1/2 ,»i (HxU n«l/2 


M x u n ) , 


(6a) 


(6b) 



while (5b) can be replaced by 


f n+1 = u n+172 + — (M u n+172 + Mu n+1//2 ) 
2 x 


(6c) 


u n+1 = v n+1 + — — (M u n+ ^ - M u n+ ^ 2 ) . 

2 y y 


(6d) 


These split forms would require additional computations, and should only be 
used to generate interpolated values at interior grid boundaries. The 
following steps illustrate one possible method of computation. 

1. Calculate v n+ ^ 72 on using (6a) 

2. Calculate u n+1 /2 on G 2 using (5a) 

3- Calculate u n+1/2 on G 1 using (6b) 

Calculate v n+1 on G 2 using (6c) 

5. Calculate u n+1 on G 1 using (5b) 

6. Calculate u n+1 on G 2 using (6d) 

Note that at each step the necessary boundary values for one grid can be 
interpolated from values at the correct level on the other grid. 

The efficacy of the alternating grid and hopscotch methods is exhibited 
in the solution of a one-dimensional model problem. The parabolic equation 

u t + (u_c)u x " ** u xx (7) 


has an exact solution 



1 2( x+1 )+(2c-1 )t 4 

u( x , t) = — (1 - tanh ). 

2 


This equation is solved on the interval [-2,2] with the exact initial value at 

t=0 and boundary values at x=-2 and x=2. A second order linearization and the 

usual central difference approximations 'are used. An overlapping set of two 

grids on the intervals [-2,. 125] and [-.125,2] is constructed. The solution 

for values of c=0. J » and p=0.05 is computed using three different methods. The 

form of the actual solution indicates that an increase in t would result in a 

translation of the graph in the positive x direction. When a numerical 

solution is computed with the Crank-Nicolson equation (2) and the values at 

x=±.125 are lagged, there is a marked deviation between the numerical and 

analytic solutions as they pass through the overlap interval. Although the 

numerical solution lags behind the actual solution, they are qualitatively 

similar with no indication of instability in the numerical solution. A 

comparison of the solutions at various times is plotted in Figure 1 . The lag 

in the numerical solution is eliminated when the alternating grid method is 

used. A careful examination of Figure 2 reveals an anomaly in the graph at 

the grid points adjacent to the interior boundary points x=±.125. This is more 

evident on the'-enlargement in Figure 3. Note that the problem occurs only at 
i' 

the points where the exceptional difference approximation is employed. The 
most accurate numerical solution for this example is calculated using the 
hopscotch algorithm. That solution appears in Figure H . In all of these 
figures, linear interpolation was used to determine solution values at grid 


boundaries. 



Methods for Conservation Laws. If the conservation equation 


u t ♦ [f(u)] x = 0 (8) 

is solved on a composite grid system, then there must be some means of 
transferring the flux f(u) from one grid to the other. There are two feasible 
alternatives. Either the solution u can be calculated by interpolation and 
then f(u) evaluated, or f(u) can be interpolated directly from the flux values 
on the other grid. The conservative difference schemes which will be 
discussed require interpolation of fluxes. However, before proceeding in that 
direction, a comparison of the two interpolation techniques will be included. 

Suppose a solution value u* at a boundary point of grid G 1 is computed by 
linear interpolation from the solution values u^ and u^ defined on grid G2. 
Then an interpolation formula of the form 

u* = au i _ 1 + gu^ a + g = 1 , 


holds, and the ^flux can be evaluated as f(u*). Now if Uq is the actual value 
of the solution at the boundary point of G-| , and hence the true flux value is 
f ( U q ) , then the interpolation procedure introduces an error as is seen in the 
following expansion. 

f(au i _ 1 +Bu i ) = f(u 0 )+f u (u 0 )(au i _ 1 +8u i -u 0 ) 



Since an expansion at the boundary point yields 


au i _ 1 + 0u i 


u o = 7 u xxo (a(x i-r x o )2 + e(x r x 0 ) 2 ), 


the leading term in the local truncation error is 


J f u (u O )u xxo (a(x i-l' x O )2 + 6(x r x 0 ) 2 ). 

Whenever the option of interpolating the fluxes is selected, then the boundary 
flux f* is calculated directly as 

f* = a r( u ± _-, ) + 

Expanding about the solution u Q , and noting the additional second order term, 


a f(u 1 _ 1 ) +6 f( Ul ) = f(u 0 ) + f u (u Q )(au i _ 1 +6u i -u 0 ) 

+ 7 f uu (u 0 )(a(u i-l" u 0 )2 + B(u r u 0 ) 2 ) 


The leading term in the local truncation error now has the form 
r' 


1 


p ^u^ u 0^ u xxo + ^uu^ u (P u xo^ a ^ x i-1 x 0 } + ^ x i x 0^ ^ 


It is clear that both procedures give an interpolation error which is 0(Ax ). 

There are many conservative finite-difference algorithms for solving 
conservation laws of the type (8). Most of the basic algorithms of practical 


interest can be written as 



u i " u i = 8( u i+1 » u ±) ~ gCuJ.uJ^). 


(9) 


For notational convenience, let 


®i+1 /2 " 8 (u i-M » u i> - 

with the implication here being that g^ + i /2 i- 3 an approximation at x i+i/ 2 * 

This notation is appropriate for the central difference approximations such as 
the Lax or Lax-Wendroff schemes. When one-sided or upwind differencing is 
used, the fractional index i+1/2 would be replaced by i or i+1 . 

Given a grid with gridpoints Xj, i=0, 1 1, the discrete conservation 

property states that 

EV = E u i + g I-1/2 " g ?/2 • 

i=l i=l 

The same result can be obtained from (8) by using numerical integration from 
x.| /2 to x j _-]/2 and the flux approximation determined by g. It is this 
derivation that will be used in the composite grid approach. 

Let G 1 and G 2 be grids defined on the intervals [a,d] and [c,b], a<c<d<b. 

y. 

The grid G 1 ha 0 points x^ i=0,1,...,I and grid spacing Ax, and G 2 has points 
y j , j=0, 1 , . . . , J and spacing Ay. The difference equations will be written in 
terms of scaled solution values v and w defined by 

v = uAx and w = uAy 


On G.j , the difference equation has the form 



and on G 2 , 



h n 

n i+1/2 


- h 


n 

i-1/2 


w n+1 

j 


„n 

W J 


w n - u n 
K j+1/2 k j-1/2 * 


There is good reason for writing the equations in this form. First of all, 
the grid spacing need not be included in the interpolation formulas, but more 
importantly, this is the required form of the difference equations when 
computing on moving grids. 

The correct interface conditions can now be derived by extending the grid 
functions to piecewise linear functions and integrating. Suppose a value k-j / 2 
is needed. Then the interval [a,b] is partitioned into two subintervals 
[ a, yi /2 ] and [y 1/2 ,b]. If y 1 /2 lies in the interval [ x i _ 1 /2 , x i+1 /2 1 , and 
h i-1/2 and h i+1 ^ 2 are known, then the value for k-jy 2 can be calculated from 
the integral property 



Assuming that and k are piecewise linear, it is easily seen that the needed 
value is the linear interpolant defined as 


k 1/ 2 - o hi- 1/2 + B h i+l/2 • 


where 


= x it M2, y 1 /? 


= yi /? _ 

x i+1/2 _x i-1/2 


a 


x i+1 /2 -x i-1 / 2 


B 



By the same argument, the interval [a,b] can be partitioned into [a,Xj_ 1 / 2 3 
and [x I _ 1/2 ,b] and the interpolation formula for the value hj_.j/ 2 * s 

h I-1/2 = a k j-1/2 + e k j+1/2 

where 

a f e = X T-1 /2~yj-1/2 # 

y j + 1/2” y j-1/2 * y j+1/2~ y j-1/2 * 


The same linear interpolation would be used on a nonuniform grid. Of 
course, the scaling factor would vary from point to point. An interpolation 
formula could also have been derived using the original equation (9), however 
the difference in grid spacing on and G 2 would have resulted in the 
appearance of a scaling factor in the interpolation formula. 

A modification of this approach can be used to develop conservative 
interface conditions for two and three-dimensional problems. The general 
two-dimensional conservation law is 

v u t + f x + 8y = °» 

1 ' 

where x and y are now the spatial variables. A difference approximation has 
the form 


.,n+1 _ _n , t.n u n , « n _ i. n / 1 a \ 

i,j = i,j h i+1/2,j h i-1/2,j k i,j+1/2 k i,j-1/2* (10) 


The grid function v is the product of the solution u and the Jacobian (or cell 
area), and the values of h and k are, up to a scalar factor, flux values in 



the direction of the curvilinear coordinate lines. Let G 1 and G2 be 
overlapping grids and suppose values of h are required along the i = 1/2 grid 
line of G 1 . For now, it is assumed that the endpoints of the grid line are on 
the boundary of the physical region. The points of G^ along the grid line are 
labeled using parameter values from any convenient parameter ization. Thus, 
let 


Pj - / 2,j * y 1 /2, j ^ * J — 0 • 1 » • • • » » 

while points of intersection of the i = 1 /2 grid line of G-| with all grid 
lines of G 2 are ordered and labeled 

q K. = ^ x i,j + 6’ y i , j +<5 ^ or ^ x i+6 , j ’ y i+6,j^’ £ =0 » 1 .**** L 

where 6 denotes a fractional index between 0 and 1 , and i, j are the indices of 
some point in G 2 . The first step in the transfer of flux values from G 2 to G-| 
is to define flux values at the points q^. If q^ lies on an i=constant grid 
line, as in the first case above, then a value hg is computed by 

interpolating the grid function k. On the other hand, if q^ lies on a 

v . 

j=constant grijd line, then hg is computed by interpolating h. Now the 
i = 1 /2 grid line divides the physical region into two parts, one covered by 
G 1 and the other covered by a subset of G 2 . If piecewise linear flux 
functions are constructed along the grid lines in each subregion and an 
integration of the flux derivatives over the complete region is performed, 
then the conservation property requires that 



y; 1 h , ♦ £jl . -5L . L £ h* * £E- . 

f-1 J 2 2 1 2 


+ 


(ID 



Here the i and n indices in (10) have been suppressed. 


The interpolation formula will be defined using a set of basis functions. 
Two additional parameter values are introduced by extrapolating from the 
parametric interval. Let q_-| = 2qQ - q.| and q^ +1 = 2q^ - q L-1 . Let ^ be the 
piecewise linear function, with knots q^, £=-1,0, , L+1, defined as 

I I , m=£ 

0 , m^£ 


where £=0,1,.. .,L, and m=-1 , 0, . . . ,L+1 . The following integrals can be easily 
computed from the parametric values of the points along the grid line. 


a £ = 


q L+1 

< 1-1 




A £’0 " 


p 1/2 

p 0 


*£ 




p j + 1 /2 
p i-1 /2 


, j — 1,2,.,.,iJ 1 


J £,J 


p J 

P J-1 /2 




These integrals are used in calculating the coefficients of the interpolation 
formulas. The formulas can now be written as 



h£ for j=0 and j=J 


h j 


L 

2 E 


i, 3 


a=o A i 


and 


= t 

1=0 




■J- hg for j=1 ,2, — , J-1 


The fact that property (11) holds is readily verified. 


11a + 


E h 

j-i 


J + 2 


J L 

-E E 

j-0 £=0 


X4- 


n 


1=0 
. M 


L i j 

E h £ 7 E 

£ j=0 


L_1 hif 

+ E h f + 

£=1 


( 12 ) 


A few remarks are sufficient to indicate how the same interpolation 
method can be employed in more general composite grid conf igurations. If 
interpolation is required on a boundary component consisting of several 
i=constant and j=constant segments, then each boundary segment could be 
treated separately with either an h value or a k value calculated from 

equations (12). However, the extrapolated parameter values would not be used 

% 

in computing the coefficients, but instead a single parameterization would be 

V 

defined for th<? entire boundary component. If the boundary component were a 
closed contour in the interior of the physical region, then the special 
boundary interpolation formulas for j=0 and j=J in (12) would be unnecessary. 

The selection of a set of piecewise linear basis functions to define the 
interpolation coefficients may be changed with only slight modification. One 
could just as easily use piecewise constant functions or use higher degree 
polynomials such as quadratics or cubics. The degree of interpolation may 
have differing effects on the numerical solution. The use of a piecewise 



constant basis may produce shock-like discontinuities, whereas a linear basis 
would tend to smear out any discontinuities in the solution. 

The conservative finite-difference scheme of MacCormack is used to solve 
the parabolic equation (7) which can be written in conservation form as 

u t + (1/2(u-c) 2 - v u x ) x = 0. 

This equation is solved on two overlapping grids on the intervals [-2,. 25] and 
[-.25,2] with the same values of 0=0.^ and p=0.05 as in the previous section. 
The numerical dissipation in the MacCormack scheme permits the use of a 
coarser grid than was used for the implicit methods. Figure 4 contains the 
solution plotted for various values of t. For this example, there was no 
noticeable difference between solutions computed with flux values interpolated 
and those computed with interpolation of solution values. 


\ 


1 



Conclusions . The accuracy of the transient solution of a hyperbolic or 
parabolic partial differential equation is dependent upon the procedures used 
to transfer information between grids in a composite grid system. The error in 
the numerical solution can be reduced by using the techniques developed here. 
While the attempt has been to construct algorithms that are easy to implement, 
the degree of difficulty would ultimately be linked to the complexity of the 


grid structure. 
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